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ABSTRACT 

We propose a new conceptual framework for the prediction of rogue waves and third-order space-time 
extremes of wind seas that relies on the Tayfun (1980) and Janssen (2009) models coupled with Adler-Taylor 
(2009) theory on the Euler characteristics of random fields. As a specific application of this framework, ex¬ 
treme statistics of the 2007 Andrea rogue wave event are examined and verified with estimates from European 
Reanalysis (ERA)-interim data. In particular, the effects of nonlinear wave-wave interactions and space-time 
variability of the wave field are considered. A refinement of Janssen’s (2003) theory suggests that in realistic 
oceanic seas characterized by short-crested multidirectional waves, homogeneous and Gaussian initial condi¬ 
tions become irrelevant as the wave field adjusts to a non-Gaussian state dominated by bound nonlinearities 
over time scales t t c ~ 0.13To/ V(T@ , where To, v and <3g denote mean wave period, spectral bandwidth and 
angular spreading of dominant waves. For the Andrea storm, ERA-interim predictions yield t c /To ~ 0(1) 
indicating that quasi-resonant interactions are negligible. Further, the mean maximum sea surface height ex¬ 
pected over the Ekofisk platform’s area is higher than that expected at a fixed point within the same area. 
However, both of these statistics underestimate the actual crest height h„i, s ~ 1.637/,. observed at a point near 
the Ekofisk site, where H s is the significant wave height. To explain the nature of such extreme, we account 
for both skewness and kurtosis effects and consider the threshold h q exceeded with probability q by the max¬ 
imum surface height of a sea state over an area in time. We find that h 0 b s nearly coincides with the threshold 
/zi /iooo ~ 1-627/j estimated at a point for a typical 3-hour sea state, suggesting that the Andrea rogue wave is 
likely to be a rare occurrence in quasi-Gaussian seas. 


1. Introduction 


Rogue waves are unusually large-amplitude surface 
waves that appear from nowwhere in the open ocean. Ev¬ 
idences that such extremes can occur in nature are pro¬ 
vided by the Draupner and Andrea events. In particular, 
the Andrea rogue wave was measured just past 00 UTC 
on November 9 2007 by a LASAR system mounted on 
the Ekofisk platform in the North Sea in a water depth of 
d = 74 m (Magnusson and Donelanj |2013| >. The Andrea 
wave has similar features of the Draupner freak wave mea¬ 
sured by Statoil at a nearby platform (d = 70 m) in January 
1995 (Haver] [2001) i. Denoting the standard deviation of 
surface elevations by o, the Andrea wave occurred during 
a sea state with significant wave height H s = 4a = 9.2 m, 
mean period Tq = 13.2 s and wavelength Lq = 220 m. The 
crest height is h = 15 m (h/H s = 1.63) and the crest-to- 
trough height H = 21.1m (H/H s = 2.3) ^Magnusson and 
|Donelan|2013| ). The Draupner wave occurred during a 5- 
hour sea state with significant wave height H s = 4ct = 11.9 
m, mean period Tq = 13.1 s and wavelength Lq = 250 m. 
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The crest height is h = 18.5 m (h/H s = 1.55) and the crest- 
to-troughheight// = 25.6 m ( H/H s = 2.15) I Haver|2004 


Magnusson and Donelan |2013| >. In the last decade, the 
properties of the Draupner and Andrea waves have been 
extensively studied iDysthe et al. 2008] Osborne |2010 


Magnusson and Donelan 2013[ |Bitner-Gregersen et ah] 


2014; Dias et al. 2015) and references therein). Several 


physical mechanisms have been proposed to explain the 
occurrence of such giant waves ( |Kharif and Pelinovsky| 
|2003| ), including the two competing hypotheses of nonlin¬ 
ear focusing due to third-order quasi-resonant wave-wave 
interactions (Janssen ( |2003| l), and purely dispersive fo¬ 
cusing of second-order waves (jFedele and Tayfun||2009| 
|Fedele|2008| i. 


For instance, recent studies propose the hypothesis that 
the Draupner wave occurred in crossing seas (see e.g. 
|Onorato et al.| ( |2010j i). These suggest that angles lying 
in the range ~ 10° — 30° between two dominant sea direc¬ 
tions are likely to lead to rogue-wave occurrences induced 
by quasi-resonant wave-wave interactions. However, |Ad-| 
|cock et af| ( |201 1} reported that the hindcast from the Euro¬ 
pean Centre for Medium-Range Weather Forecasts shows 
swell waves propagating at approximately 80° to the wind 
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sea. [Adcock et al.| ( |201 1[ > also argued that the Draupner 
wave occurred due to the crossing of two almost orthogo¬ 
nal wave groups in accord with second-order theory. This 
would explain the set-up observed under the large wave 
(W alker et al.||2004) > instead of the second-order set-down 
normally expected (Longuet-Higgins and Stewartj jl964j ). 

However, there is no evidence of significant swell com¬ 
ponents nearby the platform as clearly seen from Fig. 2 
in | Adcock et al.| ( (201 1 | ) and Fig. |T| here, which shows 
the ERA-interim wave directional spectrum at the Draup¬ 
ner site. Further, in accord with Boccotti’s (2000) quasi¬ 
determinism (QD) theory the probability that two differ¬ 
ent wave groups cross at the same point at the apex of 
their development is much smaller than the probability that 
one of the two groups focuses at the same point. One can 
also argue that reflection and diffraction from the platform 
may cause the observed set-up. However, the Draupner 
measurements were made from a bridge connecting two 
space frame structures. The structural members are rel¬ 
atively small, likely a meter in diameter. The preceding 
greatly lessens the chances for platform interference from 
spray, reflections or diffractions. Clearly, the Draupner 
wave appears fundamentally different from a typical ex¬ 
pected extreme wave because of the observed set-up of the 
mean sea level (MSL) below the large crest. However, the 
estimation of the MSL from measurements is ill-defined 
and thus may be not robust. Indeed, we note that [Walker] 
[et al.| ( |2004| estimated the mean sea level by low-pass fil¬ 
tering the measured time series of the wave surface with 
frequency cutoff f c ~ f p / 2 , where f p is the dominant fre¬ 
quency. Clearly, the mean surface and wave fluctuations 
are nonlinearly coupled and feed energetically into each 
other. As a result, the low-pass filtered mean surface is a 
mix of the two components. If the time series is not long 
enough for a statistically significant estimation of wave- 
wave interactions, the observed set-up could be the mani¬ 
festation of the large crest segment that extends above the 
adjacent lower crests. In this work, we will not dwell on 
the Draupner wave. This calls for further studies and nu¬ 
merical simulations of the sea state in which it occurred to 
clarify the physics and robustness of the observed set-up 
below the large crest. These issues are secondary and are 
beyond the scope of this paper. Indeed, we will capitalize 
on recent numerical simulations of the Andrea rogue wave 
([Bitner-Gregersen e t al.[2014[[Dias et al.|2015j ) in order to 
study the statistical properties and space-time extremes of 
the sea state in which it occurred. 

Other studies propose third-order quasi-resonant inter¬ 
actions and associated modulational instabilities as mech¬ 
anisms for rogue wave formation (e.g. Janssen] |2003j ); 
|Qsborne| ( |2010| >). Such nonlinear effects cause the statis¬ 
tics of weakly nonlinear gravity waves to significantly dif¬ 
fer from the Gaussian structure of linear seas, especially 
in long-crested, or unidirectional seas (|Janssen| (|2003|i; 
|Fedele| ( [2008| i; |Onorato et ak[ ( |2009[ ); |Shemer and Sergeeva| 


( |20091 > ; [Toffoli et aT7| ( |20 1 0| ) ; |Fedele et al.| ( |2010) ). The rel¬ 
ative importance of such nonlinearities and the increased 
occurrence of large waves can be measured by the excess 
kurtosis 

( T J 4 ) 


2*40 — 


m 


-3, 


( 1 ) 


of the mean-zero surface elevation rj, where brackets de¬ 
note statistical average. This integral statistic is defined by 


Janssen 


(2003 i as C 4 = A 40 / 3 . In general, C 4 = C 4 + C‘ 


and it comprises a dynamic component Cf due to nonlin¬ 
ear wave-wave interactions [Janssen]( j2003] >) and a bound 
contribution C? induced by the characteristic crest-trough 


asymmetry of ocean waves (Tayfun 


Lo (1990); Tayfun and Fedele (2007 

); Fedele and Tay- 

fun (2009); 

Fedele 

(2008); Janssen and BidlotJ |2009l). 


For deep-water long-crested seas with a Gaussian-shaped 
spectrum and within the framework of the higher order 
compact Zakharov (cDZ) equation (Dyachenko and Za- 
[kharov| ( |2011| >), |Fedele| ( [201 4j > showed that, correct to 
0(v 2 ) in spectral bandwidth, the dynamic excess kurto¬ 
sis monotonically increases in time toward the asymptotic 
value 


C d _ fd 

4 ,cDZ ~ ,NLS 


where 


j _ 4\/3 + ?r 2 
87T 
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4 ,NLS 
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c; 


4.NLS 


= BFI Z 


3\/3 


( 2 ) 


is the dynamic excess kurtosis of unidirectional nar¬ 
rowband waves in deep water in accord with the one¬ 
dimensional (1-D) nonlinear Schrodinger (NLS) equation 
(Mori and Janssen|( [2006[ >). Note that Eq. ([2} is also valid 
for the Dysthe (1979) equation as the associated asym¬ 
metric spectral broadening is not captured by the assumed 
symmetric Gaussian spectrum. The Benjamin-Feir index 
BFI = V2fj./v, with id denoting an integral measure of 
wave steepness defined later in Eq. ( |20l ) and the spectral 
bandwidth v is given in Eq. Clearly, C^ cDZ is smaller 
than NLS , especially as the spectral bandwidth widens. 
This is consistent with the result that in accord with cDZ 
the linear growth rate of a subharmonic perturbation re¬ 
duces with respect to the NLS counterpart for waves with 
broader spectra. Indeed, for fixed wave steepness, the 
initial-stage growth of instabilities away from a Stokes 
wave is attenuated as the spectral bandwidth increases 
|ber| ( |1978] ); |Crawford et al,| ( |1981[ )). The late-stage evolu¬ 
tion of modulation instability leads to breathers that can 
cause large waves (peregrine 19 83] |Osborne et aIT|| 2000 [ 


|Ankiewicz et a l. 2009), especially in unidirectional waves. 
Indeed, in this case energy is ’’trapped” as in a long wave¬ 
guide. For small wave steepness and negligible dissipa¬ 
tion, quasi-resonant interactions are effective in reshaping 
the wave spectrum, inducing larger breathers via nonlin¬ 
ear focusing before breaking occurs (Shemer and Sergeeva] 
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Fig. 1. Draupner storm: ERA-interim (left) directional spectrum (log scale) at the Draupner site (58.2° N, 2.5° E) at the time of maximum 
development of the storm (Jan 2nd 1995 00:00 UTC) and (top-right) wave frequency spectrum S(f)/S(f p ) and (bottom-right) angular dispersion 
g 2 D(9) = f S(f. 9)df , where cr is the standard deviation of surface elevations and f p the dominant frequency. Direction 0 = 0 means going to the 
north and 9 = 71 /2 to the east (Oceanographic convention). 


2009J; |Onorato_e^ah (2009^; [Chabchoub et al.| ( |2011 


Shemer and Alperovich (|2013|l pointed out that 


wave breaking is inevitable for ji >0.1, and breathers can 
be observed experimentally only at sufficiently small val¬ 
ues of wave steepness (~ 0.01 — 0.09) (|Chabchoub et al.| 
( |201 l|[20l2l >, see also |Shemer and Liberzon ( 2014[ l). Fur¬ 
ther, they also noted that ’breather does not breath’ and its 
amplification is smaller than that predicted by the NLS, 
in accord with the numerical studies of the Euler equa¬ 
tions ( |Slunyaev and Shrira| ( |2013| >; Slunyaev et ak|( |2013] )). 
However, typical oceanic wind seas are not 1-D but short- 
crested multidirectional wave fields. Hence, nonlinear fo¬ 
cusing due to modulational effects is diminished since en¬ 
ergy can spread directionally (Onorato et al. (2009); Tof- 
|foli et al7] ( |2010[ >). 

The sea state of the Andrea wave was short-crested (see 
Fig. [2|, and modulation instabilities may have played an 
insignificant role in the wave growth ( Alber|1978| Craw- 


bution in agreement with observations (Tayfun and Fedele 


( |2007[ >; fFedele| ( [2~008[ i; |Tayfun| ( |2008~j ); |Fedele and Tayfun 


(|2009|). This is confirmed by a recent data quality control 
and statistical analysis by |Chris tou and Ewans ( |2014| > of 
single-point measurements from fixed sensors mounted on 
offshore platforms, the majority of which were recorded in 
the North Sea. The analysis of an ensemble of 122 million 
individual waves revealed 3649 rogue events, concluding 
that rogue waves observed at a point in time are merely 
rare events induced by dispersive focusing. 

More recent studies on the statistics of extreme ocean 
waves provide both theoretical and experimental evi¬ 
dences that the expected maximum sea surface height 
over an area in time (space-time extreme) is larger than 
that expected at a fixed point (time extreme), especially 
in short-crested multidirectional seas (|Forristall (2014 


|2015|);|Fedele| ( |2012[ ); |Fedele et al.| ( |2013| i; |Barbariol et all| 


ford et al.|1981| . Further, the actual water depth to wave¬ 
length ratio, namely d/Lo ~ 0.3, suggests that waves were 
in transitional regime where modulation instabilities, if 


any, are further attenuated (see, e.g. Toffoli et al. (2009)). 
Recently, Tayfun ( 2008] ) arrived at similar conclusions 
based on the analysis of data from the North Sea. His 
results indicate that large time waves (measured at a given 
point) result from the constructive interference (focusing) 
of elementary waves with random amplitudes and phases 
enhanced by second-order non-resonant interactions. Fur¬ 
ther, the surface statistics follow the Tayfun (1980) distri- 


( 2014j >). The occurrence of an extreme in Gaussian fields 
is analogous to that of a big wave that a surfer searches 
and eventually finds (Baxevani and Rychlik i j20()6j) ). If 
the surfer searches a large area, his chances of encounter¬ 
ing the largest crest obviously increase. Indeed, a large 
number of radar signatures of likely rogue waves have 
been identified from satellite data in the context of the 


Max Wave project (Rosenthal and Lehner 2008). However, 


the relation between space-time extremes and rogue waves 
is still unclear. 

The preceding review provides the principal motiva¬ 
tion for introducing the theory of stochastic space-time 
extremes and applying it to study the space-time prop¬ 
erties of sea states in which rogue waves occur. In this 
work, we only focus on the Andrea rogue event as the 
Draupner wave requires further studies on the observed 
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Fig. 2. Andrea storm: ERA-interim (left) directional spectrum (log scale) at the Ekofisk site (56.3° N, 3° E) at about the time of the Andrea 
wave event (Nov 9th 2007, 00:00 UTC) and (top-right) wave frequency spectrum S(f)/S(f p ) and (bottom-right) angular dispersion cr 2 Z)(0) = 
f S(f , 6)df, where <7 is the standard deviation of surface elevations and f p the dominant frequency. Direction 0=0 means going to the north and 
0 = 7 t/ 2 to the east (Oceanographic convention). 


set-up ( [Walker et al.||2004] >. We aim at presenting a new 
conceptual framework for the prediction of large waves 
based on the Tayfun (1980) and Janssen (2009) models 
coupled with Adler-Taylor (2009) theory on the Euler- 
characteristics of random fields. The framework so devel¬ 
oped can be applied to study both second and third-order 
nonlinear effects and space-time properties of sea states 
where rogue waves are expected to occur. In particular, 
the new stochastic model describes the statistical structure 
of the wave surface using a Gram-Charlier type distribu¬ 
tion for crest heights formulated in ( (Tayfun and Fedele] 
2007 ) . This is able to capture the vertical crest-trough 
asymmetry induced by second-order nonlinearities mea¬ 
sured by skewness ( |Tayfun|1980j ) as well as intermittency 
effects due to third-order nonlinearities measured by kur- 
tosis (Jansseiij [2009[ >. Moreover, the statistical properties 
of local waves in space and time are related to higher order 
moments of the wave directional spectrum capitalizing on 
Adler-Taylor (2009) theory. As a result, we are able to de¬ 
fine the probability structure of the random variable rj max 
representing the maximum third-order nonlinear surface 
wave height over an area in time. In accord with the new 
model, we will show that statistical averages of the maxi¬ 
mum T] max do not explain rare occurrences of large oceanic 
waves, which instead need be interpreted using quantile- 
type approaches. In particular, we will consider the thresh¬ 
old h q exceeded by rj max with probability 0 <q<\ and the 
conditional mean h q = { rjmax |tjmax > h q \ where brackets 
denote statistical average. The statistic h q is a generaliza¬ 
tion to space-time maxima of the threshold c n exceeded 
by the crest height of the 1 jn fraction of largest waves ob¬ 


served at a point in time (Tayfun and Fedele ]2007| )). Sim¬ 
ilarly, hq generalizes the crest height average c t /„ of the 
1 /n fraction of largest waves (Tayfun and Fedele 2007 i). 

The remainder of the paper is organized as follows. 
First, the essential elements of Janssen’s (2003) formula¬ 
tion for the excess kurtosis of directional or short-crested 
seas are presented fFedele|2015[ ). This is followed by a 
review of the theory of Euler characteristics of random 
fields ( |Adler| < |T98 1 1 >), space-time extre mes (|Fedele|2012 i 
and associated stochastic wave groups (Fedele and Tayfun 
2009). We then introduce a new stochastic model for the 
prediction of space-time extremes that accounts for both 
second and third-order nonlinearities. As a specific ap¬ 
plication here, capitalizing on the ERA-interim reanalysis 
( jDee et aI7| ( |2011 ( >) and numerical simulations of the An¬ 
drea sea state (Bitner-Gregersen et al.||20l4} [Dias et al.| 


20151, the extreme statistics of the Andrea rogue-wave 


event is examined in detail. In concluding, we discuss the 
implications of these results on rogue-wave predictions. 


2. Excess kurtosis of short-crested seas 


|Fedele| ( |2015| > revisited Janssen’s (2003) formulation 
for the total excess kurtosis C 4 of weakly nonlinear gravity 
waves in deep water. This comprises a dynamic compo 


nent Q due to nonlinear wave-wave interactions (Janssen 


and Bidlot 


(2009)) and a bound contribution C 4 (Janssen 


( ]2009) >). For waves that are narrowband and characterized 
by a Gaussian type directional spectrum, Cf is expressed 
as a six-fold integral that depends on time t, BFI and the 
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parameter 


R = 


'-’e 

2v 2 ’ 


(3) 


T c = 27CV 2 ^ = —=, 

T 0 y/3R 


or 


tc 


0.13 


To VOo 


(4) 


in space and time. In the left panel of Fig. [3] the preced¬ 
ing approximation in Eq. 0 is compa red against t he the- 


which is a dimensionless measure of short-crestedness of 
dominant waves, with v and <jg denoting spectral band¬ 
width and angular spreading ( jJanssen and Bidlot] ( j2009) >; 
Mori et al. j( 201 lj >). The associated excess kurtosis growth 
rate can be solved analytically for narrowband waves 
< (Fedele| ( |20l5> , see also Appendix A). It is found that in 
the focusing regime (0 < R < 1) the dynamic excess kur¬ 
tosis initially grows attaining a maximum C 4 . max at the in¬ 
trinsic time scale 


oretical C d max for narrowband waves (Fedele (2015), see 
also Appendix A). Evidently, the latter tsTsIightly iarger 
than the maximum excess kurtosis derived by|janssen and 
|Bidlot| ( |2009[ ), who have also used the fit in Eq. ([5]t but 
with b = 1. Their maximum follows by first taking the 
limit of the excess kurtosis at large times and then solving 
the associated six-fold integral (Fedele ( 2015| >). Clearly, 
the dynamic excess kurtosis should vanish at large times. 
Janssen (personal communication, 2014) confirmed that 
Eq. ( |A3[ ) holds and provided an alternative proof that the 
large-time C d tends to zero using complex analysis. 

Further, in the focusing regime (R < 1,T C < l/>/3), 
from Eq. 0 


given by the least-squares fit 


Cj max ( R ) _ b 1 -R 
BFI 2 ~ {InfR + bRg 1 


0<1?< 1, (5) 


where Rg = and b = 2.48. Eventually the excess 
dynamic kurtosis tends monotonically to zero as energy 
spreads directionally, as in the numerical simulations of 
Annenkov and Shrira) (2009). In the defocusing regime 
(R > 1), the dynamic excess kurtosis is always negative. 
It attains a minimum at t c given by (Janssen and Bidlot 
( f2009l >) 


r d 

W.min 




4, max 


(R), 


0 </?<!. (6) 


and then tends to zero in the long time. Thus, the present 
theoretical predictions indicate a decaying trend for the 
dynamic excess kurtosis over large times. 

For time scales t > 10r c , a cold start with initial ho¬ 
mogeneous and Gaussian conditions become irrelevant as 
the wave field tends to a non-Gaussian state dominated by 
bound nonlinearities as the total kurtosis of surface ele¬ 
vations asymptotically approaches the value represented 
by the bound component (Annenkov and Shrira| (2013| 
|2014) >). In typical oceanic storms where dominant waves 
are characterized with v ~ 0.2 — 0.4 and <5g ~ 0.2 — 0.4, 
this adjustment is rapid since the time scale t c /Tg r \j 0(1) 
with Tq ~ 10 — 14 s and the dynamic kurtosis peak is neg¬ 
ligible compared to the bound counterpart. For time scales 
of the order of or less than t c , the dynamic component can 
dominate and the wave field may experience rogue wave 
behavior induced by quasi-re sonant interactions (Janssen 
(2003) )). However, one can argue that the large excess kur¬ 
tosis transient observed during the initial stage of evolu¬ 
tion is a result of the unrealistic assumption that the ini¬ 
tial wave field is homogeneous Gaussian whereas oceanic 
wave fields are usually statistically inhomogeneous both 


Q 


4,max \ ‘'w 

BFI 2 


— 1 +3t? 


(2;r) 2 1+3M 0 t 2 ' 


(7) 


Clearly, the maximum kurtosis becomes larger for longer 
time scales T c , as illustrated in the right panel of Fig. [3] 
In the defocusing regime (R> 1,T C > 1/ y/3), the dynamic 
excess kurtosis is negative and its minimum value C d min 
can be computed from Eq. (6). This result holds for deep¬ 
water waves. Drawing on Janssen and Onorato](2007j) and 
jJanssen and Bidlot) (2009) , the extension to intermediate 
waters of depth d follows by redefining the Benjamin-Feir 
Index as BFI^ = asBFI 2 , where the depth factor as de¬ 
pends upon the dimensionless depth kgd, where kg is the 
dominant wavenumber (see Appendix A). In this work we 
choose ko as the mean wavenumber k m corresponding to 
the mean spectral frequency co m = moot A w 000 > where my* 
are spectral moments (see Appendix B). In the deep-water 
limit as becomes 1, and BFIs reduces to the usual defini¬ 


tion of BFI (Janssen (20031). As the dimensionless depth 
kgd decreases, BFl^ reduces and it becomes negative for 
kgd < 1.363 and so the dynamic excess kurtosis. 

Thus, we expect that third-order quasi-resonant inter¬ 
actions do not play any role in the formation of large 
waves in realistic oceanic seas. However, the effects of 
bound nonlinearities on skewness and kurtosis should be 
accounted for in the extreme value analysis. In this re¬ 
gard, Janssen ( 2009) deeply investigated the properties of 
weakly nonlinear water waves within the Hamiltonian for¬ 
malism developed by jKrasitskii] (1994) >. Further, he de¬ 
rived general expressions for skewness and kurtosis of the 
wave surface valid for finite depths and arbitrary spectra. 
In particular, the canonical transformation for a third-order 
narrowband (NB) wavetrain can be evaluated explicitly 
and Janssen ( 2009] ) gives a simple expression for the NB 
bound kurtosis given in Eq. (A4) of Appendix A. 

Bitner-Gregersen et ak (2014) and Di as et ah) (2015[ ) 
performed Monte Carlo simulations of the hindcasted sea 
state in which the Andrea wave occurred using the irrota- 
tional and inviscid Euler equations. The estimated kurtosis 
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Fig. 3. Maximum dynamic excess kurtosis C 4 max as a function of (left) R and (right) l/r c : (bold curve) present theoretical prediction, (thin curve) 
least-squares fit from Eq. |5) (b = 2.48) and (dashed curve) Janssen-Bidlot (2009) fit (b = 1). 


is mainly due to non-resonant wave-wave interactions and 
C 4 = C\ ~ 0.033, or A 40 = A 4 6 0 = 3C 4 ~ 0.1. This is in 
fair agreement with Janssen’s NB estimate A 4( j Wfi = 0.09 
from Eq. ( |A4| i using ERA-interim data. |Bitner-Gregersen| 
|et al.| ( |20T4 i also reported a larger kurtosis (~ 0.35) sev¬ 
eral hours before the Andrea wave occurrence. The actual 
laser measurements at Ekofisk give instead oscillating val¬ 
ues within the range 0 — 0.3 indicating unstable estimates 
of kurtosis as they are based on short 20-min time series 
((Magnusson and Donela n[2013| >). 

Drawing on the ERA-interim reanalysis data, we now 
consider the hindcasted sea state during which the Andrea 
wave occurred (UTC 00 on Nov 9th, 2007). The top panel 
on the left of Fig. [4] shows the spatial distribution of sig¬ 
nificant wave height. The maximum H s is about 8.3 m, 
which is smaller than 9.2 m actually observed (Magnus- 
|son and Donelan| ( |2013| )). It is well known that ERA- 
interim underestimates peak values and prediets broader 
directional spectra because of the low space and time res¬ 
olution of the data. Indeed, wave parameters are solved 
every 6 hours and the grid cell areal size is Aq 100 2 
km 2 with 60 vertical levels (Dee et ah (2011 1 ). Neverthe¬ 
less, such predictions provide leading order estimates of 
the sea-state parameters that can be refined further in fu¬ 
ture studies, using forecast models with higher resolution 
(Ponce de Leon and Guedes Soares ( |2014) >). The top pan¬ 
els of Fig.[5]shows the Gaussian adjustment time t c /T () and 
the Janssen NB total excess kurtosis A 4 o = 3C 4 . The NB 
dynamic component A 40 = 3C 4 from Eq. ([5} and the bound 
counterpart A 40 = 3C 4 from Eq. ( |A4| > are shown in the bot¬ 
tom panels of the same figure. At the Ekofisk location the 
water depth is d = 74 m, so k^d ~ 3 and ay = 0.58. As a 


result the dynamic kurtosis is roughly half the correspond¬ 
ing value in deep waters, which is already negligible as 
the sea state is broadbanded. Clearly, the Gaussian adjust¬ 
ment time t c ~ 0(Tq) ~ 15 seconds, indicating that non¬ 
linear wave-wave interactions are negligible. Indeed, C 4 
is slightly negative, implying a defocusing wave regime 
due to the short-crestedness of the sea state whereas the 
non-zero and positive C 4 component indicates that bound 
nonlinearities are not negligible. 

Clearly, to date numerical simulations of the Euler 
equations are computationally expensive and unfeasible 
for an operational forecasting of extreme waves. One 
should consider Janssen’s general third-order expression 
for kurtosis, which is valid for any depth and arbitrary 
spectra. This is a complex formula whose implementation 
requires care and effort beyond the scope of this paper. 

In this work we propose that statistical predictions of 
extreme waves of realistic oceanic seas can be based on 
Adler-Taylor’s (2009) theory of random fields coupled 
with the Tayfun (1980) and Janssen (2009) models to ac¬ 
count for both second and third-order nonlinearities. In the 
following, we will first present the theory of space-time 
extremes for Gaussian seas ( jFedele| ( [2012[ >). Adler-Taylor 
theory is then extended to third-order non-Gaussian seas 
and apply it to study the statistical properties of the An¬ 
drea rogue wave event. 

3. Space-time (ST) extremes 

In accord with ERA-interim reanalysis, we can as¬ 
sume that sea states are both stationary in a time inter¬ 
val D < Dera = 6 hours and homogeneous over an area 
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Fig. 4. ERA-interim reanalysis at the peak of the Andrea storm. Top panels: (left) significant wave height H s — 4<7 and (right) Tayfun NB wave 
steepness jus (Eq. \2l\ ). Bottom panels: (left) wave dimension ft (Eq. ( |12| ) and (right) narrowbandedness Boccotti parameter i/r*. Dashed lines 
are H s contours. The triangle symbol indicates the Ekofisk site position. 


A < Aera = 100 2 km 2 . Then, the free surface 77 (x,t) can 
be modeled as a three-dimensional (3-D) homogeneous 
Gaussian random field over the space-time volume fl de¬ 
fined by the area A and the time interval D, and x = (x.y) 
denotes the horizontal coordinate vector. Thus, the associ¬ 
ated probability distributions at any point of the volume is 


the same and Gaussian. Drawing on Adler (1981 1 , we next 
consider the Euler characteristics (EC) of excursion sets of 
77 defined as follows. Given a threshold z, the excursion 
set Uq (z) is the part of Q. within which 77 is above z: 


Ua(z) = {(x,7) G D : 77 (x,?) > z}. 


( 8 ) 


it, plus the number of hollows inside. Further, the prob¬ 
ability of exceedance that the global maximum of 77 over 
Q, say 7 j max , exceeds z depends on the domain size and 
it is well approximated by the expected EC of the excur¬ 
sion set, provided that z is sufficiently high (jAdler|(jl981| 
|2000[ ); |Adler and Taylor| ( |2009] >). Intuitively, as z increases 
the holes and hollows in the excursion set Uq (z) disappear 
until each of its connected components includes just one 
local maximum of 77 , and EC counts the number of lo¬ 
cal maxima. For very large thresholds, EC equals 1 if the 
global maximum exceeds the threshold and 0 otherwise. 
Thus, EC of large excursion sets is a binary random vari¬ 
able with states 0 and 1 , and, for large z, 


In 1-D Gaussian processes, the EC simply counts the num¬ 
ber of z-upcrossings. Thus, © provides the generaliza¬ 
tion of this concept to to higher dimensions. Indeed, for 
two dimensional (2-D) random fields, EC counts the num¬ 
ber of connected components minus the number of holes 
of the respective excursion set. In 3-D sets instead, EC 
counts the number of connected volumetric components 
of the set, minus the number of holes that pass through 


Pr{T7 m ax >z}= Pr {EC(Ua(z)) = 1} = (EC(U a (z))), 

(9) 

where angled brackets denote expectation. This heuris¬ 
tic identity has been proved rigorously to hold up to an 
error that is in general exponentially smaller than the ex¬ 
pected EC approximation (|Adler and Taylor|(|2009|i;|Adler| 
(2000)). For 3-D random fields, which are of interest in 
oceanic applications, the probability Psj(^',A,D) that the 
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Fig. 5. ERA-interim reanalysis at the peak of the Andrea storm. Top panels: (left) Gaussian adjustment time t c /To (Eq. |JJ) and (right) total 
NB kurtosis A 40 = Ajg + Ajj). Bottom panels: (left) NB dynamic kurtosis A a() (Eq. p|) and (right) NB bound kurtosis A^q. Dashed lines are H s 
contours [m]. The triangle symbol indicates the Ekofisk site position. 


maximum surface elevation rj max over the area A and dur¬ 
ing a time interval D exceeds the threshold t;H s is given 
by ( |Adler and Taylor| ( |2009| >) 


Pst(^;A,D) = Pr {77 

max >&,} 

= (16M3^ 2 +4M 2 ^+M 1 )P R (^), 


( 10 ) 


where 


Pr(%) = Pr{/t > %H S } = exp(— 8 ^ 2 ) (11) 


is the Rayleigh exceedance probability of the crest height 
h of a time wave observed at a single point within A. Here, 
M\ and Mi are the average number of 1-D and 2-D waves 
that can occur on the edges and boundaries of the volume 
Q, and M 3 is the average number of 3-D waves that can 
occur within the volume ( |Fedele| ( |2012]> ). These all de¬ 
pend on the directional wave spectrum and are given in 
Appendix A. 

A statistical indicator of the geometry of space-time ex¬ 
tremes in the volume Q is the wave dimension J 8 defined 
by |Fedele] ( |20 1 2j > as 


4M 2 g m +2M 1 
16M 3 ^“, + 4M 2 ^ m +M! ’ 


( 12 ) 


where q m is the most probable surface elevation, or modal 
value which, according to Gumbel (1958) and Eq. 
satisfies 


P S T(§m;A,D) = (16M 3 £ 2 +4M 2 £ m +M 1 ) J P R (£ m ) = 1. 

(13) 

Drawing on jFedelej ( |2012[ ), the correct average number 
of space-time waves of dimension /j occurring within the 
space-time volume spanned by area A and time interval 
D is 


A ST (A,D) 


16M 3 ^ 2 +4M 2 ^ m +M 1 


(14) 


The parameter /j represents a scale dimension of waves, 
i.e. the relative scale of a space-time wave with respect to 
the volume’s size and 1 < j3 < 3. In particular, if wave ex¬ 
tremes are 3-D ([3 > 2) they are expected to occur within 
the volume away from the boundaries, whereas the lim¬ 
iting case of 1-D time extremes (f3 ~ 1) occur for time 
waves observed at a single point. Furthermore, |Fedele| 
(j2() 1 2)> showed that space-time extremes are larger than 
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time extremes in agreement with recent stereo measure¬ 
ments of oceanic sea states (Fedele et al7|(j2013|>, see also 

Fig. 0 . 

The bottom-left panel of Fig. [4] shows the map of the 
estimated wave dimension /3 for the North Sea area at the 
peak time of the Andrea storm for sea states of duration 
D = 3 hours over the cell area Arra- Clearly, sea states 
were short-crested and extremes are roughly 3-D, indicat¬ 
ing that the area considered is large compared to the mean 
wavelength. Thus, in accord with Boccotti’s (2000) QD 
theory, a space-time extreme most likely coincides with 
the crest of a focusing wave group that passes through the 
area as described below. 


where fl m = k m a, which is an upper bound for C 3 . From 
the linear dispersion relation k m = (0„,/g is the wavenum¬ 
ber corresponding to the mean spectral frequency (0,„ = 
mooi/mooo. 


v = xJm 0 oom 002 /mf m - 1 ( 21 ) 

is the spectral bandwidth and /n,-;* are spectral moments 
(see Appendix B). In intermediate water depths, the nar¬ 
rowband approximations lead to (Tayfun |20061 


M ~ Ms = Hmfs, 
where f$ = D\ + /A, with 


( 22 ) 


4. Extreme nonlinear wave groups 

Drawing on Boccotti’s (2000) QD theory, |Fedele and] 
|Tayfun||2009j ) showed that for second-order weakly non¬ 
linear waves the expected space-time dynamics near a 
large wave crest is that of a stochastic wave group whose 
free surface is described by 

Cc/flj =|oCt + I 0 &, (15) 


where £0 = ho/H* is the dimensionless linear crest height. 


C 1 (x,r)=¥(x,r) 

is the linear component, 

(ri{x,t)ri{x + Xj + T)) 


'F(X,T) = 


Si 


(16) 


(17) 


cos(xi)d(Oiddi 


is the space-time covariance of 77 (Bocc ottT] ( j2000 1 >) and 

C2 = / ^ Uf 2 cos (xi +Xi) + 

J ° V (18) 
A 12 cos {X\ ~X2nd(0\dQ\d(&idQ2 


is the second-order correction. Here, Sj = S((Oj,6j) 
and Xj = kj • X — (OjT, where X = (X, Y) and 
kj = (kjsin 0 j,kjcos 6 j) with kjtanh(kjd) = coj/g from 
linear dispersion, and the coefficients Ay 2 can be found 
in |Sharma and Dean] ( | 1979]> . For generic sea states, 
the largest nonlinear crest amplitude is attained at the 
focusing point (X = 0, T = 0) and given by 

? =^o + 2m?o, (19) 


where £ = h/H s is the nonlinear crest height. The Tay- 
fun wave steepness u = C 3/3 relates to the wave skewness 
C 3 of surface elevations. For oceanic applications in deep 
waters, Fedele and Tayfun] ( (20091 ) proposed the approxi¬ 
mation 

H~H a = li m ( 1-v + v 2 ), (20) 


1 

D] = -- 


4n — 1 


1 2 « 2 tanh(^ m ) — q m 

cosh(g,„) [2 + cosh( 2 < 2 ,„)] 


D 2 = 


2 sinh 3 ( q„ 


n = [1 + 2 ^r m /sinh( 2 ^ m )]/ 2 , q m = k m d and d is the wa¬ 
ter depth. The coefficients D\ and D 2 arise from the 
frequency-difference and frequency-sum terms, i.e. A 12 
and Ajj, in Eq. m as the spectral bandwidth v —» 0. A 
general expression for skewness valid for finite depths and 
arbitrary spectra is given by Janssen ( |2009( >. The NB limit 
of Janssen’s skewness yields the same Eq. ( |22| derived 
by |Tayfun] ( |2006| >. 

Further, the wave trough following the large crest oc¬ 
curs at t = T*, where T* is the abscissa of the first min¬ 


imum of the time covariance y{T) = T(X = 0,T) (Boc- 
cott il2000| . Second-order nonlinearities do not affect the 
crest-to-trough heights of large waves since wave crests 
and troughs are displaced upward equally. Thus, the max¬ 
imum second-order nonlinear crest-to-trough height ob¬ 
served at a point in time remains essentially the same 
as that of the linear group £1 and given by H/H s = 
£ 0(1 + 1 //'*), where \j/* = y(T*) is the Boccotti’s (2000) 
narrowbandedness parameter. Note that for narrowband 
waves iff* —> 1. The left panels of Fig. [4] show the maps 
of the Tayfun NB steepness u (top) from Eq. ( |22| and 
Boccotti yf* (bottom) estimated from the hindcasted ERA- 
interim sea state in which the Andrea wave occurred. 
It appears that y/* ~ 0.75 as the characteristic value of 
short-crested sea states dominated by wind waves ( |Boc-| 
c otti] ( j2000| ). At the Ekofisk location, the NB prediction 
gives /J = C 3/3 ~ 0.05 in fair agreement with the predicted 
values from numerical simulations of the Euler equations 
(Bitner-Gregersen et al.|2014 Dias et alJ2015) and actual 


TpOETli 

melanjT 


laser measurements (Magnusson and Donelan [2013| )). 

We have shown that that third-order quasi-resonant in¬ 
teractions do not play any role in the formation of large 
waves in realistic oceanic seas. However, the effects of 
third-order bound nonlinearities on the prediction of max¬ 
imum crest and wave heights must be accounted for as 
addressed in the next section. 
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5. A new stochastic third-order space-time (FST) 

model 

Drawing on|Fedele|(|2008||2012|) and|Tayfun and Fedele| 
(2007 ), we propose a new stochastic model, hereafter re¬ 
ferred to as FST, which accounts for both second and 
third-order nonlinearities in the prediction of space-time 
extremes. In particular, consider a 3-D non-Gaussian field 
over an area A for a time period of D. Clearly, the area 
cannot be too large since the wave field may not be ho¬ 
mogeneous. The duration should be short so that spectral 
changes occurring in time are not so significant and the sea 
state can be assumed as stationary. Then, the third-order 
nonlinear probability that the maximum sur¬ 

face elevation rj max over the area A and during the time in¬ 
terval D exceeds the generic threshold qH s is equal to the 
probability of exceeding the threshold £ 0 , which accounts 
for kurtosis effects only, that is 

P^-,A,D) = P ST ^o;A,D) (1 + A£ 0 2 (4 £ 0 2 - 1)). (23) 


The Gaussian probability of exceedance Pst is given in 
Eq. (10) and the amplitude (j, which accounts for both 
skewness and kurtosis effects, relates to <jjo via the Tayfun 
(1980) quadratic equation 


% =Zo + 2llZo- 
Further, the parameter 


(24) 


the wave surface relying on Janssen’s (2009) third-order 
Hamiltonian formulation, but this is beyond the scope of 
this paper. 

Given the probability structure of the wave surface de¬ 
fined by Eq. (|23l >. the nonlinear mean maximum surface or 
crest height /?fst = attained over the area A during 

a time interval D is given, according to Gumbel (1958), by 


£fst 


/?fst 

~h7 


+2p§m + 


16^ m 


7,(1 

32M 3 j m +4M 2 
1 6M 3 ^ +4 M,J m +M, 


-A G(4n)’ 
(29) 


where the most probable surface elevation value <^ m satis¬ 
fies Pp SJ ^ m ;A,D) = 1, or equivalently from Eq. (23) 

(16M 3 ^ +4 M 2 ^ m +M 1 )P R (^ m ) (1 +a£(4&- 1)) = 1, 

and 

_ 2 g m ( 8 &-l) 

iU) 1+A&(4£2-1)- 

The nonlinear mean maximum surface or crest height hj 
expected at a point during the time interval D follows from 
Eq. (29) by setting Mi = M 3 = 0 and M\ = Ad, namely 

§r = h T / H s = +2/x^ + + ’ ( 3 °) 

l0Sm—AG(5 m ) 


A = A 40 + 2 A 22 + A 04 


(25) where, now, satisfies 


is a measure of third-order nonlinearities and it is a func¬ 
tion of the fourth order cumulants X nm of the wave surface 
r ] and its Hilbert transform rj (Tayfun and Fedele 20071 

Krn = (t]"f) m ) /o n+m + (—l) m/ 2 (n— l)(m — 1 ), (26) 

where a is the standard deviation of the wave surface and 
n + m = 4. Drawing on Mori and Janssen)( [20061 >, we as¬ 
sume the following relations between cumulants 


A D /^ m )(l+A^(4<^— 1} ) = 1 . 

Here, Ad = D/T denotes the number of wave occur¬ 
ring during D and T is the mean up-crossing period (see 
Appendix B). The linear mean counterpart follows from 
Eq. (30) by setting = 0 and A = 0. 

Drawing on the Boccotti (2000) distribution for wave 
heights, the third-order nonlinear mean maximum wave 
height expected at a point is given by 


A22 = A40/3, 


A 04 = A 40 , 


(27) 


fl T =*T A V2(l + r), (31) 


and Eq. (|25]> simplifies to 


A — A appr — 


8A40 

~~3~~ 


where is the nonlinear mean maximum crest height 
that accounts for kurtosis effects only. This follows from 
(28) Eq. (30) by setting ju =0, namely 


which will be used in this work. Then, Eq. (23) reduces to 
a modified Edgeworth-Rayleigh (MER) distribution ( |Mori| 
and Janssen 2006). For realistic oceanic seas the kur¬ 
tosis A 40 ~ A*,, is only due to bound nonlinearities. To 
date, the validity of the cumulant relations in Eq. (27) has 
been proven to hold for second-order NB waves only |Tay-| 
|fun and L o 1990). Further studies are desirable to deter¬ 
mine the general expressions of fourth order cumulants of 


h W 


T /Hs ~ U+ 


(32) 


where 


A D ^(^ m )(l+A^(4^_i)) = L 


For narrowband waves, i/P = 1 and Eq. a reduces to 
the MER model proposed by jMori and Janssen (2006). 
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Clearly, this overestimates wave heights in short-crested 
or multidirectional seas as, in general, yf* < 1 . 

When the lateral dimension t = VA is much larger than 
the average wavelength Lq, the maximum surface height 
occurs most likely within the area of interest and not on 
its boundaries. In this case, on average the number of 3- 
D waves is much larger than the numbers of 1-D and 2-D 
waves that can occur on the boundaries, i.e. M 2 2> M 2 and 
M\. Keeping only the leading term in Eq. ( |29[ > yields the 
following asymptotic value of the expected surface wave 
height maximum over large areas 

s(3D) 

*=FST 

where 


h i3D) 

rl FST 

H, 


— 


Ye (1 +4/i<^n 


16 £ m -f-AG(£ n 


(33) 


1.6 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



16 M 3 CPr(U) (1 + A^(4^ - 1)) = 1. 

For 2-D waves only (Mi = M 3 = 0 in Eq. (29) ) 

y e (l+4ju£ m ) 


.(2D) 

^s? ) = %=^ + 2 




16 §m-f -AG(§m)’ 
(34) 


where 


4 M 2 ^ m P R (U) (1 + A^(4^ - 1)) = 1. 


If one ignores kurtosis effects (A = 0), Eqs. ([33) and (34) 
reduce to the 2-D and 3-D analogues of the Piterbarg- 


Tayfun (PT) model formulated by Socquet-Juglard et al. 
( 2005) , and hereafter referred to as 2D-PT and 3D- 
PT respectively (see also [Piterbarg] (T995 ) and Forristall| 
(20TT) ). 

In offshore applications, the interest is in the expected 
wave maxima over small areas such as those covered by 
oil rigs, i.e. I < L {] . In this range, |Fedele et al.](2013) have 
shown that the boundary corrections accounted by both 
terms M\ and M 2 are important for a correct estimation of 
Euler characteristics and expected maxima (see also For- 
|ristall|(20T5)). Hence, they cannot be ignored as has been 


assumed by Romolo and Arena (2015), since maximum 
surface heights expected over small areas are underesti¬ 
mated. This point is elaborated further in the next sec¬ 
tion and demonstrated explicitly by way of the results dis¬ 
played in Fig. [7] 

As an application, consider now the hindcasted sea state 
during which the Andrea wave occurred (wave steepness 
M=C 3 /3 0.05 and excess kurtosis A 40 = 3 C 4 ~ 0.1 from 

Dias et al.| (20l5 >). The top panel of Fig. [ 6 ] displays the 


predictions from Eq. (30) for the third-order mean max¬ 
imum surface or crest height hj expected at a point as a 
function of the excess kurtosis A 40 for a typical D = 3- 

hour sea state. From Eq. 


32 




the mean height h\ that ac¬ 
counts for kurtosis effects only is also shown. Comparing 


Fig. 6 . Andrea wave: (top panel) nonlinear mean maximum surface 
or crest height hj expected at a point as a function of excess kurtosis 
A 40 , which accounts for both second-order (skewness) and third-order 
(kurtosis) nonlinearites. The prediction h^' of average heights that 
accounts for third-order kurtosis effects only is also shown. (Bottom 
panel) nonlinear mean maximum wave height expected at a point (Hj 
) as a function of excess kurtosis A 40 . Dashed lines denote the actual 
Andrea crest and wave height values (h os j, = 1.63//, and fi a b s = 2.3 H s ). 


the two predictions, it is clear that second-order nonlinear¬ 
ities cannot be neglected as they yield a substantial 15% 
increase in crest height in comparison to the modest 5% 
increase due to kurtosis (~ 0 . 1 ) from the linear estimates 
(A 40 = 0). Note that even for unrealistic values of kurtosis 
(~ 0 . 8 ) the predictions are well below the observed actual 
Andrea crest height h„b s = 1.63 H s (dashed line). Further, 
the bottom panel of Fig. [ 6 ] displays the mean maximum 
wave height Hj expected at a point from Eq. (31 1 . This 
is smaller than the observed actual Andrea wave height 
Hobs = 2.3 H s . 

In summary, averages of maxima do not explain the ac¬ 
tual large Andrea time crest and wave heights observed at 
the measurement site. This point will be discussed and 
addressed later in section [7] 


6 . ST analysis of the Andrea wave 

In this section we will study the space-time properties 
of the Andrea storm in comparison to stereo wave mea¬ 
surements at the Acqua Alta site of a short-crested sea 
state dominated by bora winds (experiment 2, see |Fedele| 
|et al.|(20l3) ). For example, Fig.[7]shows the prediction of 
space-time extremes at the Acqua Alta site and from the 
hindcasted ERA-interm Andrea sea state during which the 
rogue wave occurred. In particular, the left panel of the fig¬ 
ure shows the Acqua Alta observed ratio of space-time and 
time maximum surface heights over a period of D ~ 0.5 
hours as a function of the lateral dimension £ normalized to 
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Fig. 7. Left panel, Acqua Alta stereo measurements (Fedele et al.|j201~3) ): comparison between the observed ratio of space-time and time 
maximum surface heights (hollow squares) and theoretical averages /trsT/ ^T as a function of the lateral side length i normalized to the average 
wavelength Lq. Fedele space-time (FST) model (solid line) and 2D and 3D Piterbarg-Tayfun (PT) models nearly the same as 2D- and 3D-FST 
(dashed lines). Right panel, Andrea ERA-interim predictions: expected ratio of mean areal and point maximum surface heights estimated at the 
Ekofisk location (FST, black line; 3D-FST, dashed line). For comparison the estimated ratios for Acqua Alta are also shown (FST. blue solid line; 
3D-FST 2D-PT, blue dashed line). 


the average wavelength Lq ~ 19 m (hollow squares, max¬ 
imum length i max ~ 24 m). In the same panel, we com¬ 
pare the theoretical FST ratio /(i stAt, which accounts 
for boundary corrections associated with M\ and M 2 in Eq. 
( |29| , and the asymptotic 3D-FST from Eq. ( |33j ), which is 
valid over large areas and ignores boundary effects. The 
2D-FST model from Eq. ( |34| ) is also shown. They both 
nearly coincide with 2D-PT and 3D-PT as the measured 
kurtosis at Acqua Alta was almost Gaussian. As one can 
see, the boundary contributions cannot be neglected over 
areas with lateral dimension comparable to or smaller than 
the typical wavelength. Indeed, both 2D- and 3D-FST un¬ 
derestimate the observed ratios when l < Lq. Note that for 
£ ~ Lo the observed space-time wave surface maximum 
/ifst is 1.4 times larger than the time maximum hj at a 
single point. 

A similar trend is also observed for the expected space- 
time extremes of the Andrea sea state for the same time 
period of D = 0.5 hours. Specifically, the right panel of 
Fig. [7] displays the mean maximum surface height ratios 
/ifst/^t (FST, black curve) as afunction o££/Lq estimated 
at the Ekofisk location, where the ERA-interim predicts a 
mean wavelength Lq ~ 150 m. Note that the FST ratios 
for Acqua Alta and Andrea (blue and black solid curves) 
are close to each other for £ < Lq. These results are very 
encouraging as the observed sea state at Acqua Alta was 
in deep waters (measured H s = 1.09 m and d/ Lq = 1.25, 
with d = 16 m) whereas the Andrea sea state was in in¬ 
termediate waters (ERA-interim estimates H s = 8.3 m and 
d/Lo = 0.49, with d = 74 m). Although ERA-interim un¬ 


derestimates the actual value of H s at Ekofisk (= 9.2 m), 
it appears that the maximum surface height ratio /ifst/^t 
is slightly sensitive to the significant wave height level and 
just depends on average spectral properties of the sea state. 
Further studies are desirable to investigate possible sta¬ 
tistical similarities and universal laws for space-time ex¬ 
tremes in wind sea states, but this is beyond the scope of 
this paper. 

Clearly, space-time extremes cannot explain the actual 
large Andrea time crest height observed at the measure¬ 
ment point. The space-time analysis of the Andrea storm 
simply predicts that the maximum surface wave height 
over the platform footprint area is ~ 20 % higher than the 
average maximum surface height that is expected at a fixed 
point within the same area (t/L q ~ 0.3, see right panel 
of Fig. [7]). Indeed, in relatively short-crested directional 
seas such as those observed around the Ekofisk area, it 
is very unlikely that the observed crest actually coincides 
with the largest crest of a group of waves propagating in 
space-time. In contrast, in accord with Boccotti’s (2000) 
QD theory, it is most likely that the sea surface was in fact 
much higher somewhere near the measurement point. 


7. Threshold h q exceeded with probability q by the 
maximum surface height 

Consider the average duration of D = 3 hours for the 
Andrea sea state in which the rogue wave was observed, 
considering that the ERA time resolution is of Dera = 6 
hours. Monte Carlo simulations yield wave steepness fi ~ 
0.05 and excess kurtosis A 40 ~ 0.1 (Bitner-Gregersen et al. 
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|2014[|Dias et al.|2015| l. It follows that the mean maximum 
nonlinear wave surface height hj expected at a point is 
~ 1.2 H s (see Fig. (§)■ This is lower than the actual value 
hobs ~ 1-63 H s observed. Statistically, this means that in an 
ensemble of M realizations of 3-hour sea states, each of 
which has similar statistical structure to the Andrea wave 
field, all the maximum surface heights observed at a point 
will exceed 1.2 H s . Clearly, the maximum surface height 
at a point can reach or exceed the actual observed value 
1.63 H s only in few realizations out of the ensemble of M 
sea states. 

To characterize such rare occurrences in weakly nonlin¬ 
ear third-order random seas one can consider the thresh¬ 
old h q = i % q H s exceeded with probability 0<<?<1 by the 
maximum surface height rj max over an area A during a sea 
state of duration D. This satisfies (see Appendix C) 


P^ q -A,D)=q, 


(35) 


where £, q = £o, 9 +2 J u£d 0 > £ 0 , q follows from Eq. (C5i 




and the nonlinear probability of exceedance P IST is given 
in Eq. <1231. Further 
(t]max|t?max > h q ) is given by 


the conditional mean h q = 


hq — h q -p {Xi.q A PX2,q TA ^3 q A ^4 q)H s / (36) 


where the coefficients Xj,q depend on q as well as on the 
spectral parameters Mj and q q . For Gaussian waves, the 

linear threshold h q L> and conditional mean ]} q ] follow by 
setting both fi = 0 and A = 0 in Eqs. ( [35] ) and ( [36] respec¬ 
tively. For second-order waves the associated threshold 

hq 1 * and mean !n q 1 follow by setting A = 0 only. As q be¬ 
comes smaller the threshold h q tends to coincide with the 


mean h q since terms Xj^/Q i n Eq. (36i tend to zero, and 
similarly for linear and second-order thresholds. 

In this work, we consider the time thresholds h^ q , h^ q 

and hj q at a point and conditional averages hj q , li-y q and 
hj q . These can be computed using Mj = M 3 = 0 and 
Mi = A)) in Eqs. ([35} and ( [36] . Here, Ad denotes the 
number of zero-crossing waves occurring during the time 
interval D. Further, note that /it , q =\ coincides with the 
mean maximum third-order surface height hj and simi¬ 
larly for linear and second-order conditional averages, that 

is hj j = hj and h^l = h^\ 

Drawing on the Boccotti (2000) distribution for wave 
heights, the third-order nonlinear threshold exceeded with 
probability q by the maximum wave height at a point is 
given by 


H^ = h^y/2{l + r), 

and the conditional mean 


(37) 


H T . q = h^y/2(l + r)- 


(38) 


Here, the threshold lij q accounts for kurtosis effects only 
and it follows from Eq. ([35] with jj. = 0. Similarly, the 




conditional mean h T follows from Eq. (36 


sian waves, the linear counterparts are given by 


For Gaus- 


H. 


(G 


T,9 


= 4 ^ 2 ( 1 +r) 


(39) 


and 

H^ q =h^2{\ + r)- ( 40 ) 

Note that Hj q= 1 coincides with the mean maximum non¬ 
linear wave height Hj and the linear conditional mean 

ijL _ tjL 

n j j — n r p. 

The statistical interpretation of the threshold h q and 
conditional mean h q is as follows. Consider M realiza¬ 
tions of a stationary and homogeneous sea state of du¬ 
ration D. On this basis, there would be M samples, say 
(t?max, • tlmax) of the maximum surface height r] max ob¬ 
served within an area A during the time interval D. Then, 
the threshold h q is exceeded in only qM realizations of the 
ensemble of M sea states, and h q is the average of the qM 
largest values in the sample. Similar interpretation holds 
for the threshold H q and conditional mean H q for the max¬ 
imum wave height. 

For the time statistics at a point, hj q coincides with the 
standard threshold c„ exceeded by the crest heights of the 
1 /« fraction of largest waves if we set n = Ad/ q, where Ad 
is the average number of zero-crossing waves expected in 
a D-hour sea state. Similarly, lr\\ q is the crest height aver¬ 
age c\/ n of the q/N D fraction of largest time waves (see, 
for example, |Tayfun and Fedele| ( |2007} ). Clearly, the pre¬ 
vious interpretation does not hold for space-time or mul¬ 
tidimensional random fields. Indeed, the notion of zero- 
upcrossings is generalized to that of Euler Characteristics 
of excursion sets and wave counting is not as obvious as 
in one-dimensional random processes as there is no clear 
definition of a space-time wave ( |Fedele|2012} . 

Consider now the hindcasted ERA-interim Andrea sea 
state during which the rogue wave occurred. Fig. [ 8 ] dis¬ 
plays the third-order threshold hj q from Eq. (36} for the 
maximum crest height at a point as a function of the ex¬ 
ceedance probability q for a 3-hour sea state (thin black 
line). The linear hj (thin red line) and second-order 

threshold (thin blue line) are also shown in the same 
figure. The corresponding conditional averages are also 
displayed (bold lines). If we account for both second 
and third-order nonlinearities, the actual crest height point 
measurement observed (/?„/„ ~ 1.63 H s ) nearly coincides 
with/i x . 1/1000 ~ 1.62 H s . Thus,onlyin 1 sea state out of the 
ensemble the maximum surface height exceeds 1.62 H s . 
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Fig. 8 . Andrea wave: third-order nonlinear threshold hj ^ for the 
maximum crest height at a point as a function of the exceedance prob¬ 
ability q for a 3-hour sea state (thin black line). The linear threshold 
/*Tr/ (^in red line) and the second-order threshold (thin blue line) 

are also shown. The conditional averages h q , and are also 
displayed (bold lines). The horizontal dashed line indicates the actual 
Andrea wave crest value. Note that coincides with the mean 

maximum surface height hj and similarly for linear and second-order 
conditional averages, i.e hj j = hj and h^\ = The wave steep¬ 
ness /.i ~ 0.05 and the excess kurtosis A 40 ~ 0.1 |Bitner-Gregersen et al.| 
[20l4l|Dias et al.|2015) . 

For second-order nonlinearities only, h 0 b s is slightly be¬ 
low ^ 1/10000 ~ l-68//y (thin blue line). Further, h 0 t, s ex¬ 
ceeds the linear ^J/kjoooo ~ 1-52 H s and it nearly coin¬ 
cides with hj j/j 0Q0 00Q making the Andrea wave event an 
extremely rare occurrence in Gaussian seas. 

The predictions for the maximum wave height at a point 
are shown in Fig. & Note that the actual wave height ob¬ 
served at Andrea (H 0 b s ~ 2.3 H s ) exceeds the third-order 
threshold //t,i/io = 2.13 H s and it nearly coincides with 

the linear Hj^/jgg = 2.26 H s . This indicates that the An¬ 
drea wave most likely occurred in the configuration of a 
maximum crest and the associated wave height was not 
the largest and not as extreme as the crest height. Indeed, 
according to Boccotti’s (2000) QD theory, extreme waves 
with the largest wave height have almost symmetrical 
crests and troughs, but the Andrea wave had a strong crest- 
trough asymmetry (Magnusson and Donelan]( |2013) >). 

Note that h T i/^ooo i s the same as the threshold cpooo.ooo 
exceeded by the crest of one wave in a sample of n = 
No/q = 10 6 waves, where q = 10~~’ and Ad ~ 10 3 (mean 
period ~ 12 s). Thus, in ideal stationary conditions that 
last forever the Andrea crest height h 0 b s is exceeded in 
average every R = 0.5 years. However, oceanic seas are 
non-stationarity and a long-term analysis of storms is re¬ 
quired to estimate the probability of occurrence P s of a 
D-hour sea state as that similar to the Andrea field. Then, 


the probability that the threshold h q is exceeded during the 
lifetime span L of an offshore structure is P = P s q, which 
is an important design parameter. The associated return 
period R follows from the relation P = exp (-L/R), which 
assumes a Poisson statistics for storm occurrences ( jBoc-[ 
c otti|2000] ). Such an analysis requires predicting extremes 
of non-stationary wave fields following |Haver| (j2002j) (see 
also ( |Boccotti|2000||Fedele|2012] >), but this is beyond the 
scope of this paper. 


In summary, second-order skewness nonlinearities of 
crest heights are dominant as also indicated by extensive 


statistical analyses of oceanic wave measurements (Tay- 

fun and Fedele 

2007, Tayfun 2008; Fedelej|2008 

Fedele 

and Tayfun 2009). Further, in accord with Janssen 

2003), 


bound kurtosis effects cannot be neglected and must be ac¬ 
counted for in order to obtain more accurate estimates of 
both crest and wave height extremes. 



Fig. 9. Andrea wave: (top) linear and nonlinear thresholds H^ and 
H„ (thin lines) of the maximum wave height at a point as a function of 
the exceedance probability q for a 3-hour sea state. The conditional av¬ 
erages H'* ] and H q (bold lines) are also shown. The horizontal dashed 
line indicates the actual Andrea wave height value. Note that Hj q ] 
coincides with the mean maximum wave height fly and similarly the 
linear conditional average fly | = Hj. The excess kurtosis A 40 ~ 0.1 
(Bitner-Gregersen et al.|2014||Dias et al.|20151. 


8 . Concluding remarks 

Current freak wave warning systems rely on Janssen’s 
(2003) theory for the kurtosis of surface elevations, a key 
result with significant implications to the understanding 
of the role of nonlinear wave interactions in rogue wave 
formation ( |Janssen and Bidlot|2009| . The present results 
suggest that in typical oceanic fields third-order quasi¬ 
resonant interactions do not appear to play a significant 
role in the wave growth. Fedele’s (2015) refinement of 
Janssen’s theory shows that the large excess kurtosis tran¬ 
sient observed during the initial stage of wave evolution 
is a result of the unrealistic assumption that the initial 
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wave field is homogeneous and Gaussian. Oceanic wave 
fields are typically non homogeneous both in space and 
time, and initial conditions become irrelevant as the wave 
field tends to a more realistic non-Gaussian nonlinear state 
where statistics are affected mainly by bound harmonics 
( (Annenkov and Shrira| ( |2013| |2014| >). In this regime, sta¬ 
tistical predictions of extreme waves can be based on the 
Tayfun (1980) and Janssen (2009) models to account for 
both second and third-order nonlinearities. In particular, 
skewness effects on crest heights are dominant and bound 
kurtosis must be accounted for to attain more accurate pre¬ 
dictions of extremes. 

Our statistical analysis of the Andrea storm reveals that 
space-time extremes cannot explain the actual large An¬ 
drea time crest and wave height observed at the measure¬ 
ment point. The ST analysis simply predicts that the mean 
maximum sea surface height expected over the Ekofisk 
platform’s area is higher than the mean maximum height 
expected at a fixed point within the same area, irrespective 
of the significant wave height level. The actual maximum 
occurs at some point within the area or its boundaries, and 
the likelihood of that point’s coincidence with the point 
where measurements are done is essentially nil. 

Further, both of the time and space-time averages of 
maxima underestimate the actual Andrea crest and wave 
heights observed. If one accounts for both second and 
third-order nonlinearities, the actual crest height nearly co¬ 
incides with the threshold /? i /iooo exceeded with probabil¬ 
ity 1 /1000 by the maximum crest height at a point in a typ¬ 
ical 3-hour sea state. This suggests that the Andrea rogue 
wave is likely to be a rare occurrence in quasi-Gaussian 
seas. Instead, the actual wave height is not so extreme 
as the crest height as it nearly coincides with the thresh¬ 
old //wjoo exceeded with probability 1/100 by the max¬ 
imum wave height at a point. The wave height was not 
the largest as wave measurements indicate a strong crest- 
trough asymmetry of the Andrea wave, and according to 
Boccotti’s (2000) QD theory this most likely occurred in 
the configuration of a maximum crest. 

Finally, the present conceptual framework for wave ex¬ 
tremes can be applied to higher order resolution wave fore¬ 
casts and operational models for more accurate predictions 
of rogue waves. 


analysis of Acqua Alta measurements and Guillermo Gal- 
lego for his support with LaTeX. FF acknowledges partial 
support from NSF grant CMMI-1068624. 


APPENDIX A 


Dynamic and Bound Kurtosis 

For narrowband waves in deep waters, the evolution of 
the dynamic excess kurtosis from initial Gaussian condi¬ 
tions is given by ( jFedele|f2015| >) 

C 4 = BFI 2 J(x.R) (Al) 


where 

J(v,R) =2Im 


1 


da. 


0 V1 - 2/a + 3a 2 \/l + 2iRa + 3 R 2 a 2 

(A2) 

The maximum is attained at T = t c (see Eq. (|4|) and given 

by 


Ci mm (R)=BFI z J„(R ), 


(A3) 


where 


Jp(R) = J 


= j( 




:R 


= Im 


, 1 
V3R 


0 x/1 -2/a + 3a 2 \/l +2iRa + 3R 2 a 2 


da. 


and Im(cz) denotes the imaginary part of a. 

Drawing on Janssen and Onorato ( |2007| ) and Janssen 
and B idiot (2009), the dynamic kurtosis in intermediate 
waters of depth d is simply computed by replacing the 
Benjamin-Feir Index with 


BFlj=BFI 2 a s 


where the depth factor 

a - ( c s \ 2 S x,, l 

S \c 0 ) kpCOpCO^ 

depends on the mean frequency (Op 

(Op = \Zgk 0 T 0 , D 0 = tanh (k 0 d) 
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corresponding to the mean wavenumber ko via the linear 
dispersion relation, the group velocity 

/ 1 f 2kod 1 ©0 

c, = © 0 = - C o|l+ sinh(2M) |, co=-, 

where a>p is the first derivative of the angular frequency 
with respect to the wavenumber ko, and the second deriva¬ 
tive 


®0 = ~g 


{Dp — kod (l — Dq)} 2 +4 (kpd) 2 Pf ) (l — Dp) 
4©o kpDp 
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Further, the nonlinear interaction coefficient 

■>4 inn2 in i I />*- _ /">\2 


_ 9Dq — IOD 5 + 9 1 

nl ml kod 


1 + 


( 2 c g -c 0 / 2 )- 
r 2 _ r 2 


where c s = \fgd is the phase velocity in shallow waters. 

The NB bound kurtosis C 4 is given by 
in the form 


Janssen 


(2009) 


C 4 — A 40/3 — 8 /jg (/3 + y+2(a + A)"). 
where the wave steepness jlQ = k^o. 


a = 


3 ~ T o 2 
4T () 3 ’ 


13 = 


3 8 + (l -7^) 


2\3 


64 


r r 6 
J 0 


r= ~2 a ~’ 


and the wave-induced mean sea level variation 


respectively. They all depend on the mean period T, mean 
wavelengths L x and Ly in a and y directions: 


T=2nJ — , L x = L y = 2nJ m ' m 


mo 02 


m 200 


mo 20 


(A4) 


and 


Here, 


a„, = 1 - al - a 2 , ~ a 2 ,+ 2a xt a yt a xy . 


m uk=[\ k i x kif k s(f,e)dfde 

are the moments of the directional spectrum S(f, 6 ) and 


1 / — 


mioi 


: • OCy, = 


m 011 


mn 0 


yjm 2 oomoo2 v /mo20"Jo02 ' v / '”200'mo20 


A 


1 


c 


2 

s 


1 — T 1 

A_4 

T 0 


1 

x 


APPENDIX C 


In deep waters, C\ = 6 /i ( 2 . Note that Eq. ( |A4| ) is not valid 
for small water depths as second and third-order terms of 
the associated Stokes expansion can be larger than the lin¬ 
ear counterpart (see Eq. (A18) in Janssen| f2009] >) . This 
implies the constraints aj,i m < 1 and < 1 , which 

are not violated for the Andrea sea state. Note that fur¬ 
ther restrictions apply as occurrences of spurious crests on 
the troughs of large, relatively steep second-order Stokes 
waves are anomalous and not an inherent characteristic of 
real waves ( |Tayfun|2013| >. 


Threshold h q exceeded with probability q by the 
maximum surface height 

For 3-D random fields, which are of interest in oceanic 
applications, the nonlinear probability P^j(^\A,D) that 
the maximum surface elevation q max over the area A and 
during a time interval D exceeds the generic threshold qH s 
is equal to the probability of exceeding the threshold £ 0 , 
that is 

P$&A,D) =Pst(^o;A,D) (1 +A£o 2 (4£o 2 - 1)), (Cl) 


APPENDIX B 

Space-Time Statistical Parameters 

For space-time extremes, the coefficients in Eq. ( p~4| ) are 
given by (Baxevani and Rychlik| ( |2006] >; [Fedele] < |2012[ >) 

„ , D 4 ty 

Mi=2n fTTZ a X yt ’ 


where the Gaussian probability of exceedance (see Eq. [TO] 
and ( jAdler| 198T[ |Adler and Taylor|2009[ >) 

Ps T (£o;A,D)=Pr{77 

max > SoH s } 

= (16M 3 ^ 0 2 +4M 2 ^o+A7i)exp(-8^ 0 2 ), 

(C2) 


the parameter A is given in Eq. 28 and amplitudes § and 
£0 are related by the Tayfun (1980) quadratic equation 


M 2 =\/2jt (= = ^JT— 


at, 


D L r x 4 l y r -— 

- = = \/l - a 2 , + = = - a 2 

T Ly V Ly Ly V ^ 


Ml =N D +N x +Ny, 


where 


D 

Nd = =, 


A 7 

N x = =, 

T 

‘-‘x 


N v = =4 


t, =4, + 2 J u^o 2 - 


(C3) 


The nonlinear threshold h q = E, q H s exceeded with proba¬ 
bility q by T] max is given by 


are the average number of waves occurring during the time 
interval D and along the x and y sides of length 4 and £ y , 


pW(S q -,A,D)=q, 

or equivalently 

P ST (^, q -,A,D) (1 + A^^o 2 , - 1 )) = q, 

where from Eq. ( |C3| ) 

= ^0,q+2/1^0,q. 


(C4) 


(C5) 


(C 6 ) 
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The mean h q of j] max conditioned on rj ma x > h q then fol¬ 
lows by definition from 


h q /H s 


f^P ( ^(4lA,D)d4 

fgpM(5;A,D)d? ’ 


(C7) 


where the nonlinear pdf 




d_P^_ 

dE, 


(C8) 


Integrating by parts yields 


h q /H s = 




,(nl) 


-JZPg’MdS 
& 


where 

poo 

Xi 4=1 PsT&dE, 

■'SO ,q 

= pi, 9 exp(-8^ 0 2 i9 ) + + M 3 )Erfc(2v / 2£ 0 , 9 ), 

poo 

X2,q = / 4 yP ST {y)dy 

4 SO ,q 

= P2,q ex P( — 8 ^o 9 ) + ^^M 2 Erfc(2v/2£ 0 , ? ), 

poo 

X-i , q = / y 2 (4y 2 -1) P ST ( y)dy 

J %0,q 

= p 3 , 9 exp(-8^ 0 2 9 ) - M ‘ 51 ^ M3 v / 27rErfc(2v / 2^ 0 ^), 

poo 

X 4 ,q = / 4y 3 (4y 2 - l) P ST {y)dy 

o _ 

= p4 !9 exp(-8§ 0 2 ? ) + ^-M 2 Erfc(2V2^o, 9 ), 

Erie (a) is the complementary error function and 


and it follows from Eq. ( |C4| ) that 

hq/Hs = Zq+ l 

4 ■'S? 


Next, note that 


PpST&dt, = I P ST (&) (1 + A<^o 2 (4^T - 1)) dE ,, 




(C9) 

where E,q (E, ) is a function of the nonlinear E, given by 
Eq. (C3i. A change of variables via E,o{E,) = y and 
dE, = ffifidy = (1+4 fly)dy simplify the preceding expres¬ 
sion to 


rP s ( f (§ )d% = r Pst (y)F(y)dy, (CIO) 

J £q ■'lo ,q 

where F(y) = (1 + Ay 2 (4y 2 — 1)) (1 +4 fly). Further, §o i9 
follows from Eq. ( |C5| i and it relates to via Eq. ( |C6| >. 
Therefore, 


hq/Hs = Zq + ~ [ Psr(y)F(y)dy. 

C 1 ,q 

This expression is written in the form 


hq — hq + (xi,q T PX2,q “E A ^ 3 ^ + [J.A.%4 q'j H s /q, (Cl 1) 


Pi ,q =4 (M 2 +4M 3 ^o.q), 

P2, 9 =7 (44i +2M 3 +4M 2 ^o, 9 + 16M 3 ^^ q 


P2,q —^0 ,q 


3M 3 — M\ M\ + M 3 e 2 ,3 „4 


64 
3M 2 


■ 3 7 £o, , q +^ 2 ^o ,o +4M 3 i ® 0 


_M 3 

P4 ' 9_ l6 + 




M2 f 3 4 . 

■ t5o+ 


(Mi +2M 3 )^o i9 + 4M 2 ^o ? + 16M 3 ^o ? . 
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